The stationary density matrix of a pumped polariton system 
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The density matrix, p, of a model polariton system is obtained numerically from a master equa- 
tion which takes account of pumping and losses. In the stationary limit, the coherences between 
eigenstates of the Hamiltonian are three orders of magnitude smaller than the occupations, meaning 
that the stationary density matrix is approximately diagonal in the energy representation. A weakly 
distorted grand canonical Gibbs distribution fits well the occupations. 
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The possibility of Bose-Einstein condensation (BEC) 
of excitonic polaritons in microcavities has raised big ex- 
pectations recently [U |3] . Due to the relatively small 
lifetimes of polaritons (of the order of picoseconds) and 
still smaller rates for phonon relaxation in the lower po- 
lariton branch [1], a very important question should be 
answered concerning whether the observed magnitudes 
come from a system in thermal equilibrium or have a 
dynamical nature. 

In paper [5], we show that in a decaying system the 
enhancement of ground-state occupations can be un- 
derstood in terms of the combined effects of polariton- 
polariton scattering and photon emission, even if phonon 
relaxation does not act. This suggests that the results 
of the experiment reported in |T] could be ascribed to a 
dynamical effect and not necessarily to BEC in the po- 
lariton system. 

On the other hand, in the continuously pumped sys- 
tem, where a stationary state is reached when pump and 
losses are equilibrated, it was undoubtedly demonstrated 
that phonon relaxation is not effective in the lower polari- 
ton states [6]. Thus, the question arises about what kind 
of stationary state are we reaching in the experiments 
reported in [2]. 

In the present paper, we are aimed at giving a par- 
tial answer to the latter question. We will assume that 
there are not thermalization mechanisms, and will com- 
pute the density matrix arising from a purely dynamical 
equation. We use the same model of polariton system, 
previously studied in Refs. |5l |7j, with a finite num- 
ber of single-particle states for electrons and holes and 
a single photon mode in the microcavity. A term ac- 
counting for pumping is added to the master equation 
for the density matrix. This master equation is numer- 
ically solved in order to find the stationary density ma- 
trix. The main result of the paper is the following: in the 
stationary limit the density matrix is approximately di- 
agonal in the energy representation and, thus, describes a 
kind of quasiequilibrium which can be fitted to a weakly 
distorted grand canonical Gibbs distribution. The dis- 
torted Gibbs distribution can be thought of as coming 
from the maximization of the entropy with an additional 



constraint in phase space commuting with the Hamilto- 
nian and fulfilling the requirement of additivity. This 
idea was presented in Rcf. [8; to describe quasi station- 
ary nonequilibrium states and share similarities with the 
results of Hamiltonian dynamics simulations [9J . 

For completeness, we first recall the main features of 
the model polariton system. The Hamiltonian describing 
the system is the following: 

i 

ijkl ijkl 

- (3^{ij\\kl) elhXhjek + {Egap + A) a^a 

ijkl 

+ g^^a)hje, + ae\h\^ . (1) 

i 

We include 10 single-electron and 10 single- hole levels 
in Eq. ([ij (the first three two-dimensional harmonic- 
oscillator shells) . The single-particle spectrum for elec- 
trons and holes is assumed fiat: 

T'>'^=Egap. Tj'^^^O. (2) 

This means that the dot confinement energy, hcoo, is as- 
sumed much smaller than the effective band gap, Egap. 
Our model describes a relatively small quantum dot 
strongly interacting with the lowest photon mode of a 
thin microcavity. (3 is the characteristic Coulomb en- 
ergy, P = e'^ /{Anelosc) = / (ine) ^ mujQ / h, where e is 
the medium dielectric constant, and lose = ^/^-/(towo) - 
the oscillator length. We will take the value, (3 =2 meV. 
{ij\\kl) are matrix elements of the Coulomb interaction 
between harmonic oscillator states. The parameter A 
gives the detuning of the photon energy with respect to 
the (bare) pair energy, equal to Egap, and g =3 mcV 
is the photon-matter coupling strength. Notice that the 
photon couples the electron state i to the hole state i, 
which differs from the latter only in the sign of the angu- 
lar momentum. This means that electron-hole pairs are 
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FIG. 1: (a) The lowest states in sectors with < Npoi < 
10. (b) The matrix elements |(F|a|/)p, smeared with a 
Lorentzian of width F = 0.1 meV, vs the energy difference 
E{I) - E{F) - Egap. The states \I) correspond to Npoi = 2, 
whereas |F) is the ground state in the Npoi = 1 sector 



created or annihilated in states with zero angular momen- 
tum. As a result, the Hamiltonian, Eq. Q, commutes 
with the total angular momentum of the system, Ltotai- 
In what follows, we consider only Ltotai = states. 

The Hamiltonian, Eq. ([l]), preserves the polariton 
number, 



Npoi = a^a + y2{h\h-, + e\e,)/2. 



(3) 



We diagonalize the Hamiltonian in a basis constructed 
from Slater determinants for electrons and holes and Fock 
states of photons. For a given polariton number, Npoi, 
the wave functions arc looked for as linear combinations: 



1^) = X! ^S,,SH,n\^ej Sh, 



(4) 



where and Sh are Slater determinants with the same 
number of particles, Npairs, and n = Npoi—Npairs- When 
Npoi — there is only one state, the vacuum. When 
Npoi = 1 there are 17 states with Ltotai = 0. One of them 
is the state with one photon (no pairs), and the remaining 
16 states correspond to matter excitations (no photons), 
that is, all possible combinations of one electron and one 
hole states with total angular momentum equal to zero. 



As Npoi increases, the number of eigenstates of H rises, 
reaching around 18000 for Npoi = 10. We use Lanczos 
algorithms jim to obtain the energies and wavefunctions 
of the lowest 20 states, which are used to write down the 
master equation for the density matrix. This number of 
states, 20, is chosen on the basis of two reasons. First, 
we want to keep the number of matrix elements, ppi, 
between reasonable limits. And, second, in a sector with 
given Npoi we shall include all of the states, |/), with 
significant matrix elements, {F\a\I), between |/) and the 
ground state in the sector with Npoi — 1- They are very 
important in the dynamics, as it will be shown below. 

We draw in Fig. [ija) the set of states used in this work 
for a detuning A = — 3 meV. Notice the approximate 
linear dependence of the ground-state energy on Npoi, 
and the energy gap from the ground to the first excited 
state in each sector, roughly proportional to ^/Np^i g. 

A sample of the matrix elements |(F|a|7)p is repre- 
sented in Fig. [ijb). Matrix elements are smeared out 
with a Lorentzian of width F = 0.1 meV. The states |/) 
correspond to Npoi — 2, whereas \F) is the ground state 
in the Npoi = 1 sector. The lower (LP) and upper polari- 
ton (UP) branches are clearly distinguished. Notice that 
the latter are included among the 20 states which partic- 
ipate in the dynamics. This is important because they 
substantially contribute to the population of the lowest 
state in each sector with given Npoi- 

Now we come to the central point of the paper, the 
master equation for the density matrix and its stationary 
solution. The master equation is written as [Tl] [T2]: 



dp 
dt 



P 



- — [H, p] + — {2apa) — a) ap — pa) a) 
^{2a\jpaij - aijaljp - paija\j). (5) 



The parameter k accounts for photon losses through the 
cavity mirrors {hn « Egap/Q, where Q is the cavity qual- 
ity factor). In our calculations, we take k — 0.1 ps^^. 
Notice that k << g, thus our model system works under 
the strong light-matter coupling regime. Other sources 
of losses such as, for example, spontaneous exciton decay 
are much less important and will not be considered. On 
the other hand, the parameter P is a pumping rate. In 
our modeling, incoherent pumping is supposed to come 
from highly excited excitonic states, which decay towards 
the lower polariton states through phonon emission. This 
is a kind of polariton reservoir. We will use a sort of 
homogeneous pumping, with equal probabilities for all 
states. To this end, we introduce lowering and rising op- 
erators, cr/j|/) = I J), crljl J) = |/). As we are employing 
a finite number of states (20) in each sector with given 
Npoi > 1, total pumping probabilities are finite. Notice 
that we are not including relaxation acting at the "hor- 
izontal" level in Fig. [T]^a), i.e. causing the excited po- 
lariton states to decay towards the lowest states with the 
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FIG. 2: (a) The occupations in each polariton sector. The 
detuning parameter is A = —3 meV, and the pumping rate 
is P — 0.01 ps^^. Ground-state occupations are indicated by 
squares, (b) The best fit to the actual occupations with a 
grand canonical Gibbs distribution, (c) Best fit with a dis- 
torted grand canonical Gibbs distribution. 



same Npoi. Thus, if the lowest polariton states become 
occupied in our simulation it is the result of dynamics and 
not of relaxation. The absence of phonon thermalization 
is also the reason why Ltotai = states are decoupled 
from other states with Ltotai 7^ 0. 

In the stationary limit, the r.h.s. of Eq. ^ is equal 
to zero. This set of homogeneous linear equations can be 
shown to be linearly dependent. Indeed, it can be easily 
verified that 



d ^ 



(6) 



which corresponds to the conservation of probability, 
^iPii = 1- We replace the equation for pn by the con- 
straint pii = 1 in order to obtain an inhomogeneous 
system: 



J,K 



MpijK Pjk = Bpi, 



(7) 



where Mn jj = 1, the rest of the matrix elements, 
AIpijK, are obtained from Eq. (Sj), and the components 
of the vector Bpj are zero, with the exception of Bn = 1. 
Only the occupations, p//, and the "horizontal" coher- 
ences, pfi, where F and / are states in the same Npoi 
sector, acquire nontrivial values as a result of solving Eq. 
([t]). The other coherences are equal to zero. 

We show in Fig. |2ja) the resulting occupations for 
A = —3 meV and P — 0.01 ps""'^. The mean number 
of polaritons is (Npoi) — 4.77. In this case, the system 
is under the "polariton laser" regime [6 . That is, above 
threshold {{Npoi) ~ 1) and below saturation due to Fermi 
statistics {{Npoi) much greater than 10 in the present 
model, which is called "photon laser" regime in Ref. i6j). 

The first point in the figure, enclosed by a square, cor- 
responds to vacuum's occupation. The next 17 points 
corresponds to states with Npoi — 1- Then, there are 20 
states corresponding to Npoi = 2, etc. In order to fa- 
cilitate lecture of the figure, we have indicated the first 
state in each sector by enclosing it with a square. We 
may notice that the ground-state occupations in sectors 
with Npoi < {Npoi) are notably enhanced as compared 
with the rest of the states in the same sector. This 
is mainly the result of a transfer of population from 
"upper-polariton" states, following from the matrix ele- 
ments illustrated in Fig. [ijb). Consequently, the "upper- 
polariton" occupations are depressed. 

The computed horizontal coherences are three orders 
of magnitude smaller than the occupations. For A = —3 
meV and P — 0.01 ps~^, we obtained J2i<j \Pij\ — '''•^x 
10^^, which should be compared with pn = 1. This 
means that the stationary p is approximately diagonal in 
the energy representation. This is the main result of the 
paper. We verified that the above statement holds for 
any values of the system's parameters. 

Once shown that H and p approximately commute, 
we can try to answer the question about whether the 
computed p can be obtained from a maximum entropy 
principle. In a first attempt, we fit the occupations to 
a grand canonical Gibbs distribution, p ~ exp(— (i? — 
PiNpoi)/Ti). For the situation represented in Fig. |2ja), 
we get the effective parameters, Ti = 28.05 meV, and 
Pi = Egap — 11.7 meV. The r.m.s. deviation is 0.003, one 
fourth of the maximum values of p. Qualitatively, how- 
ever, it gives a poor description of the actual occupations, 
as it can be seen in Fig. |2]^b). 

A second possibility for the fit is motivated by the re- 
sults of paper where the stationary p of a nonequilib- 
rium system is obtained from the maximization of the en- 
tropy with an additional constraint in phase space. Then, 
we should find the extremum of the functional: 
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P (ps-^) 


/i2 - Egap (meV) 


Ta (meV) 


02 (meV ^) 


0.001 


-2.60 


13.92 


1.4 xlO~^ 


0.01 


-0.51 


8.37 


1.5 xlO~^ 


0.1 


-0.10 


4.00 


6.9 xlO"'' 



TABLE I: Effective parameters corresponding to the fit given 
by Eq. (flO|. A = -3 meV. 



S = -^pnln{pji) + aiC^pii - 1) 
I I 

+ 02 PiiE{I) -{E))+ a^iQip) - (Q)) 
I 

+ a4{J2Pii^poiiI)-i^poi))^ (8) 

where ai, . . . , 04 are Lagrange multipliers, and Q{p) is 
the constraint. It was argued in Ref. [5] that the require- 
ments of additivity and commutativity with the Hamilto- 
nian fix the form of the constraint, Q{p) = J2i P^IE{I), 
where g is a kind of Tsallis index, and F commutes with 
p. We choose F = H — EgapNpoi- The equation for p, 

= - Mpii) - 1 + ai + a2E{I) + a^Npoiil) 
+ a^qp]-\E{I) ~ EgapNpoiil)). (9) 

under the assumption that the distortion is weak, 9 ~ 1, 
can be iterativcly solved. Taking the first iteration, we 
get: 

The parameter 012 = a^q^q—l) is expected to be small. In 
our model, fitting the computed occupations to Eq. ( 10 1, 
we obtain T2 — 8.42 meV, p2 = Egap — 0.53 meV, and 
a2 — 0.0015 mcV^^. The r.m.s. deviation is 9 x 10^^, 
three times smaller than in the previous case. Quali- 
tatively, the obtained distribution, represented in Fig. 
|2jc), is much more closer to the actual occupations. The 
analogs of lower and upper polariton states show the 
highest dispersion. 

In Table |T] we give the effective parameters for pump- 
ing powers in the interval from 0.001 to 0.1 ps~^. In these 
computations, we neglect coherences and extend Eqs. ([t]) 
for the occupations up to a maximum Npoi — 40. We no- 
tice that the dependence of the chemical potential on P 
is qualitatively the same as observed in experiments jT3] . 
The effective temperature decreases with P, but for still 
higher values of the pumping starts increasing, as in the 
experiments. We stress, however, that the experimen- 
tal fits cover a range in the emission angle, whereas our 
L = distribution is qualitatively related to fc = states, 
responsible for the emission along the normal direction. 



In conclusion, we have computed the stationary den- 
sity matrix of a pumped polariton system, which comes 
from a dynamical master equation without thermaliza- 
tion mechanisms. The density matrix is shown to be ap- 
proximately diagonal in the energy representation. Thus, 
a kind of quasiequilibrium is established. Although 
our model describes a quantum dot strongly interact- 
ing with the lowest photon mode of a thin microcav- 
ity, a quasiequilibrium of dynamical origin could also be 
present in the experiments described in Refs. [U [SJ [3]. 
We should notice in this respect that the similarity be- 
tween lasing and a second-order phase transition was un- 
derlined long ago in the context of laser physics (see for 
example Ref. [121 Chapter 11, and references therein). 
The accumulated evidences on Bose-Einstein condensa- 
tion of polaritons could be a present-day manifestation 
of this old idea. 

A dynamical framework could be also the basis for the 
computation of the photoluminescence spectra, second- 
order coherence functions, etc in the present model |14j . 
Experimental facts such as the low threshold for polari- 
ton lasing, increase of linewidth with pumping power, etc 
are nicely (qualitatively) reproduced. In this context, our 
work for a multiexcitonic quantum dot O [71 [M] explores 
an intermediate (in the number of states) region between 
the single-level dot [TT| and the infinite system (quantum 
well, see for example Refs. [T5]). where the system 
is simple enough to be studied by exact diagonalization 
methods, but complex enough to exhibit many of the 
properties of the infinite system. 

We showed that the computed density matrix can be 
reasonably fitted to a weakly distorted grand canonical 
Gibbs distribution. The multi-polariton analogs of UP 
and LP branches are the states worst fitted. Let us stress 
that the same strategy of obtaining p from entropy max- 
imization with an additional constraint in phase space 
can be apphed [T7], with success, to the description of ex- 
periments on metaequilibrium states in electron plasma 
columns [18'. 
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